

%% data, voting- and campaign-stage estimates, expected surpluses and winning probabilities (and gradients)

load('voting_stage.mat')
load('campaign_stage.mat')
load('coalition_stage_p1.mat') % for "short" replication
% coalition_stage_p1 % for "long" replication, uncomment this line


%% "entry" parameters: parsimonious

% no-candidate payoff parameters
X_PVEM = [ones(N,1),zeros(N,1),PW0_PVEM];
X_PRI = [zeros(N,1),ones(N,1),PW0_PRI];
optCFS = optimset('Display','off','HessFcn',@(thetaNC,lambda)HessLL_cfs(thetaNC,lambda,X_PVEM,X_PRI,ES0,ES1,ES2));
[theta_noCovs,~,exit_cfs_noCovs] = knitro_nlp(@(thetaNC)LL_cfs(thetaNC,X_PVEM,X_PRI,ES0,ES1,ES2,[N0,N1,N2]),...
    randn(size(X_PVEM,2),1),[],[],[],[],[],[],[],[],optCFS,'knitro.opt');

% robust inference (accounting for voting- and campaign-stage estimation uncertainty)
% gradient of joint moment conditions
G_all = grad_cfs_moments(theta_noCovs,X_PVEM,X_PRI,ES0,ES1,ES2,DES0,DES1,DES2,DPW0_PVEM,DPW0_PRI);
G_all = [G_all;... % gradient of coalition-stage MLE FOCs
    zeros(size(Ghat,1),length(theta_noCovs)),Ghat,zeros(size(Ghat,1),length(B2S0));... % gradient of campaign- and voting-first-stage MCs
    zeros(size(X2S0,2),length(theta_noCovs)+length(gamma)+length(B)),X2S0'*X2S0]; % gradient of voting-second-stage MCs
% joint covariance matrix
% coalition-stage score
P1 = exp(X_PVEM * theta_noCovs + ES1 - ES0);
P2 = exp(X_PRI * theta_noCovs + ES2 - ES0);
P0 = 1 + P1 + P2;
P1 = P1 ./ P0;
P2 = P2 ./ P0;
L = (([zeros(N0,1);ones(N1,1);zeros(N2,1)] - P1) .* X_PVEM + ...
        ([zeros(N0+N1,1);ones(N2,1)] - P2) .* X_PRI)';
L = - L';
% campaign-stage MCs
Zc = (A(1:N,:) ./ sum(A(1:N,:),2)) * (Z .* (C - br0));
% voting-first-stage MCs
Zxi = (A(1:N,:) ./ sum(A(1:N,:),2)) * (Z0 .* xi0);
% voting-second-stage MCs
Zxi2S = [zeros(N0,size(X2S0,2));(repmat(A(N0+1:N,N0+1:N),1,2) / 2) * (X2S0 .* xi2S0)];
Omega_all = [L' * L,L' * Zc,L' * Zxi,L' * Zxi2S; ...
    Zc' * L,Zc' * Zc,Zc' * Zxi,Zc' * Zxi2S; ...
    Zxi' * L,Zxi' * Zc,Zxi' * Zxi,Zxi' * Zxi2S; ...
    Zxi2S' * L,Zxi2S' * Zc,Zxi2S' * Zxi,Zxi2S' * Zxi2S];
V_all_noCovs = inv(G_all' * (Omega_all \ G_all));
seTheta_noCovs = sqrt(diag(V_all_noCovs));
seTheta_noCovs = seTheta_noCovs(1:length(theta_noCovs));
rTheta_noCovs = [theta_noCovs,seTheta_noCovs,2*(1-normcdf(abs(theta_noCovs./seTheta_noCovs)))];

disp(' ')
disp('for Table 6 (column (I)):')
disp(print_results_theta(round(rTheta_noCovs,3),{},{'constant'}))


%% "entry" parameters: rich

% no-candidate payoff parameters
X_PVEM = [dChar.region==1,dChar.region==2,dChar.region==3,dChar.region==4,dChar.region==5,D,zeros(N,5+KD),PW0_PVEM];
X_PRI = [zeros(N,5+KD),dChar.region==1,dChar.region==2,dChar.region==3,dChar.region==4,dChar.region==5,D,PW0_PRI];
optCFS = optimset('Display','off','HessFcn',@(thetaNC,lambda)HessLL_cfs(thetaNC,lambda,X_PVEM,X_PRI,ES0,ES1,ES2));
[theta,~,exit_cfs] = knitro_nlp(@(thetaNC)LL_cfs(thetaNC,X_PVEM,X_PRI,ES0,ES1,ES2,[N0,N1,N2]),...
    randn(size(X_PVEM,2),1),[],[],[],[],[],[],[],[],optCFS,'knitro.opt');

% robust inference (accounting for voting- and campaign-stage estimation uncertainty)
% gradient of moment conditions
G_all = grad_cfs_moments(theta,X_PVEM,X_PRI,ES0,ES1,ES2,DES0,DES1,DES2,DPW0_PVEM,DPW0_PRI);
G_all = [G_all;... % gradient of coalition-stage MLE FOCs
    zeros(size(Ghat,1),length(theta)),Ghat,zeros(size(Ghat,1),length(B2S0));... % gradient of campaign- and voting-first-stage MCs
    zeros(size(X2S0,2),length(theta)+length(gamma)+length(B)),X2S0'*X2S0]; % gradient of voting-second-stage MCs
% joint covariance matrix
% coalition-stage score
P1 = exp(X_PVEM * theta + ES1 - ES0);
P2 = exp(X_PRI * theta + ES2 - ES0);
P0 = 1 + P1 + P2;
P1 = P1 ./ P0;
P2 = P2 ./ P0;
L = (([zeros(N0,1);ones(N1,1);zeros(N2,1)] - P1) .* X_PVEM + ...
        ([zeros(N0+N1,1);ones(N2,1)] - P2) .* X_PRI)';
L = - L';
% campaign-stage MCs
Zc = (A(1:N,:) ./ sum(A(1:N,:),2)) * (Z .* (C - br0));
% voting-first-stage MCs
Zxi = (A(1:N,:) ./ sum(A(1:N,:),2)) * (Z0 .* xi0);
% voting-second-stage MCs
Zxi2S = [zeros(N0,size(X2S0,2));(repmat(A(N0+1:N,N0+1:N),1,2) / 2) * (X2S0 .* xi2S0)];
Omega_all = [L' * L,L' * Zc,L' * Zxi,L' * Zxi2S; ...
    Zc' * L,Zc' * Zc,Zc' * Zxi,Zc' * Zxi2S; ...
    Zxi' * L,Zxi' * Zc,Zxi' * Zxi,Zxi' * Zxi2S; ...
    Zxi2S' * L,Zxi2S' * Zc,Zxi2S' * Zxi,Zxi2S' * Zxi2S];
V_all = inv(G_all' * (Omega_all \ G_all));
seTheta = sqrt(diag(V_all));
seTheta = seTheta(1:length(theta));
rTheta = [theta,seTheta,2*(1-normcdf(abs(theta./seTheta)))];

disp(' ')
disp('for Table 6 (column (II)):')
disp(print_results_theta(round(rTheta,3),dnames,{'region1','region2','region3','region4','region5'}))


%%

clearvars -except optCFS theta_noCovs exit_cfs_noCovs ...
    V_all_noCovs seTheta_noCovs rTheta_noCovs ...
    X_PVEM X_PRI theta exit_cfs V_all seTheta rTheta

save('coalition_stage.mat')



